QSGW calculation of the work functions of Al(lll), Al(100), and Al(llO) surfaces 
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Modifications to the quasiparticle self-consistent GW (QSGW) method needed to correctly de- 
scribe metal/ vacuum interfaces and other systems having extended regions with small electron 
density are identified and implemented. The method's accuracy is investigated by calculating work 
functions for the Al(lll), Al(100), and Al(110) surfaces. We find that the results for work function 
do not depend on the DFT functional employed to calculate the starting Hamiltonian and that 
QSGW yield results in quantitative agreement with data from ultrahigh vacuum experiments. 
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I. INTRODUCTION 

The work function is not only a most important quan- 
tity characterizing the surface of a metal, it also di- 
rectly affect surface phenomena like growth rate, the form 
of crystallites, sintering, catalytic behavior, adsorption, 
chemical reactions, surface segregation, and formation of 
grain boundaries. In addition, the work function largely 
determines rates of electron surface emission and is there- 
fore of applied interest for optimizing thermionic emit- 
tersi, where a low work function is sought, and pulsed- 
power components, e.g, in transmission lines, where a 
high work function is required. An ability to make accu- 
rate theoretical calculations of work functions is therefore 
of great interest and importance for the development of 
new high-performing materials. 

Presently, most calculations of work functions are 
made within density functional theory (DFT^i 3 - using 
either the local density approximation (LDA) or the gen- 
eralized gradient approximation (GGA). Unfortunately, 
the accuracy of the DFT-based methods in calculation of 
the work function is not always satisfactory and results 
depend on the DFT functional used (see, e.g., recent re- 
view of available theoretical and experimental values of 
work function for a number of metals in Ref. 0) ■ For in- 
stance, the work functions of Al surfaces calculated by us- 
ing the GGA with Pcrdcw-Burke-Ernzcrhof (GGA/BPE) 
functional^ are approximately 0.2 eV lower then cor- 
responding LDA values calculated with Ceperley- Alder 
(LDA/CA) functional^ and approximately 0.3 eV lower 
then corresponding LDA values calculated with Barth- 
Hedin (LDA/BH) functional 7 - (see Table I). This exam- 
ple highlights the need for methods capable of reaching 
0.1 eV (or better) accuracy for theoretical prediction of 
the work function of metals. 

The GW approximation of Hcdin 8 is a well estab- 
lished method which yields highly accurate quasiparti- 
cle (QP) energies for bulk materials 9,10 . GW calcula- 
tions are usually performed in a non-self-consistent man- 
ner by using the LDA Green's function, G, and screened 
Coulomb interaction, W (so-called GqWq method). The 



results obtained by the Go Wo method depend on the 
quality of underlying DFT wave-functions and eigenval- 
ues and the agreement with experimental data worsen 
if the DFT description is not sufficiently accurate. For 
example, Go Wo fails to describe the band gap of NiO 11 . 
Recently, Faleev, van Schilfgaarde and Kotan i 10 ' 11 de- 
veloped the so-called Quasiparticle Self-consistent GW 
(QSGW) method which is independent of DFT and 
demonstrated that results for QP energy levels of bulk 
materials obtained by the QSGW are in better agree- 
ment with experiment then results obtained by the stan- 
dard GqWq method. In contrast to the GoWo method, 
QSGW describes correctly also strongly correlated ma- 
terials like NiO or MnG±i. 

Although calculations of surface QP energies were per- 
formed already more then two decades agoi£, GW calcu- 
lations for surfaces and other non-bulk systems remains 
rare due to the demanding computational requirements. 
Most commonly, the GW method has been applied to 
study an image potential and corresponding image states 
of insulating^ 3 , and metallic*^— surfaces and clusters^ 7 -. 
Recently, the GW method was used to investigate the 
image potential-induced renormalization of the molec- 
ular electronic levels for molecules adsorbed on metal 
surface o 18 ' 19 and thin insulator films 2 ^. Note that the 
image potential cannot be obtained by the LDA or GGA 
approaches since they do not include non-local polariza- 
tion effects that are present in GW theory through the 
W operator. Another application of the GW method 
to non-bulk materials that is becoming an active area 
of research is QP calculations of transport properties of 
nanoscale systems^t-S 3 .. 

To the best of our knowledge, two GW studies of the 
work functions of metals have been published to date. 
Morris et al2i calculated the work functions for Al(lll), 
Al(100), and Al(110) surfaces using the G W method 
and making a jellium approximation. They also included 
vertex corrections to the self-energy and W (evaluated 
on a homogeneous electron gas level) that resulted in 
a significant, over 1 eV, underestimation of calculated 
work functions, attributed to inherent self-interaction er- 
ror. Heinrichsmeier et al 2 ^ proposed a new non-local 
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parametrization of the exchange-correlation functional 
derived from the Go Wo calculations for jellium surfaces 
(the Vz C (GW) method), and applied the method to real 
(111) and (100) surfaces of Al and Pt. A conclusion 
of both studies was that the values of the Al(lll) and 
Al(100) work functions obtained by the GoWo method 
are significantly worse then corresponding LDA values as 
compared to the experimental data (see Table I). Since 
both these GW calculations involve the jellium approxi- 
mation, the question remains how well the fully atomistic 
GW method could describe the work functions of metals 
as compared to the DFT results and experiment. 

In this paper, we present results of work function cal- 
culations for the Al(lll), Al(100), and Al(110) surfaces 
evaluated within fully atomistic GqWq and QSGW ap- 
proaches. As we show below, the values of the work 
function obtained by these two methods only insignifi- 
cantly deviate from each other (within 0.02 eV), do not 
depend on the initial DFT functional, and are in excellent 
agreement with experimental data. On the other hand, 
the DFT approaches predict different results for the work 
function depending on the functional used. The paper is 
organized into two main sections describing the methods 
and the results followed by a concise summary. 

II. METHOD AND COMPUTATIONAL 
DETAILS 

Let us first briefly describe the computational ap- 
proach. QSGW is a method to determine nonlocal (but 
static and Hcrmitian) optimum one-particle Hamiltonian 
H° in a self-consistent wa y 10 ' 11 ' 26 . First, starting with a 
trial Hamiltonian H° (usually, the LDA Hamiltonian is 
used as the first-iteration H°) the self-energy £(w) is cal- 
culated in the GW approximation. The static self-energy 
is defined in the basis of the eigenfunctions V'kn(r) of the 
Hamiltonian H° as follow a 10 ' 11 ' 26 

= i?e(^ k „|[£(e k „) + £(ew)]/2|Vw) . (1) 

where k is the wave vector, n is the band index, £ kra de- 
note eigenvalues of H , and Re means to take the Hermi- 
tian part. Next, the Hamiltonian H° is updated for each 
iteration using S„ n / instead of the usual LDA exchange- 
correlation potential 

H° = ~ + +v H +Y J \^n) SL'Wwl , (2) 

knn' 

where V ext is the external (nuclei) potential and V H is 
the Hartree potential. This procedure is iterated until 
self-consistency is reached. Here, self-consistency is de- 
fined as when £^ n / generated by H° is identical (within 
a small tolerance) to the Ej^, that enters into H°. It 
has been showni - that this procedure in an approxi- 
mate way minimizes the difference between the full non- 
local, non-static and non-Hermitian GW Hamiltonian 

H ( U} ) = + yeXt + V " + E H and Hamiltonian H° 



(that is why it is called 'optimum'). Note that H(ui) 
is a functional of H° because both V H and E(w) cal- 
culated in the GW approximation depend on eigenfunc- 
tions generated by H°. Hence, the iteration procedure 
described above self-consistently determines both H(u>) 
and the corresponding optimum H°. The Go Wo method 
is simply the first iteration of described above cycle, self- 
energy is obtained from Eq. (fT|) using LDA wave func- 
tion and energies, and first-iteration Go Wo Hamiltonian 
is constructed by Eq. @. 

In the present work, we used the experimental lat- 
tice constant of Al at zero temperature, 4.025 A2£. The 
Al(lll), Al(100), and Al(110) surfaces were modeled by 
a (1 x 1) surface unit cell in xy-directions and a peri- 
odic combination of Na Al layers and Ny vacuum lay- 
ers in z-direction. The vacuum layer was of the same 
width as the Al layer, and contains so-called floating ba- 
sis orbitals 2 ^ placed instead of the atomic muffin-tin or- 
bitals. Na ranged from 4 to 12 and Ny ranged from 6 
to 10 were used to analyze the convergence of the results 
with respect to these parameters. The work function, 
is defined as the difference between the electrostatic 
potential at a point far from the surface and the Fermi 
energy; $ = V es (oo) — sp. In our calculations V es (oo) 
was estimated as electrostatic potential in the middle of 
the vacuum slab. 

It is known that GW calculations of band gaps in 
semiconductor thin film s 29 ' 30 and molecular chains 3 -^ per- 
formed in repeated-cell geometries converge slowly with 
vacuum thickness due to the long-range nature of the 
non-local screened Coulomb interaction. Rozzi et al^ in- 
vestigated this problem and developed a Coulomb cut-off 
scheme to eliminate the long-ranged slab-slab interaction. 
They found that the introduced Coulomb cutoff parame- 
ter mostly affects delocalized unoccupied states. On the 
other hand, both quantities that enter the expression for 
the work function: the Fermi energy and electrostatic po- 
tential are affected only by occupied states that are local- 
ized inside the metal slab and therefore only weekly de- 
pend on the thickness of the vacuum slab. Consequently, 
the work function converges quickly with the number of 
vacuum layers: we found that Ny=6 is sufficient to de- 
termine the work function to 0.005 eV accuracy. 

Both relaxed and unrelaxed internal atomic coordi- 
nates were utilized to model the Al surfaces. The relaxed 
position of Al layers with Na > 10 were taken from GGA 
calculations of Da Silva 32 . In particular, for Al(lll) sur- 
face the values of interlayer expansion/contraction were 
+1.15%, -0.05%, +0.46%, +0.21%, -0.05% (first num- 
ber is for the surface layer, last number for fifth layer), 
for Al(100) surface +1.59%, +0.44%, -0.02%, -0.68%, - 
0.56%, and for AlfllO) surface -7.18%, +3.87%, -2.12%, 
+2.04%, +0.82% [32[. 

The QSGW method is implemented as an exten- 
sion of the all-electron full-potential linear muffin-tin or- 
bital (LMTO) program suite. The diagram of the self- 
energy and density/potential self-consistency cycles that 
includes the LMTO and GW parts of the code are shown 
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on the Figure 2 of Ref. [26J . The description of the basis 
sets and other details of the LMTO and QSGW imple- 
mentations can be found in Refs. (2(| and [28T ]. 

The surface Brillouin zone (BZ) integration in the 
LMTO part of the code were performed with (22 x 22) 
and (24 x 24) Monkhorst-Pack meshes^ (k LMTO -mesh) . 
The GW self-energy was calculated with (6x6), (8x8), 
and (10 x 10) meshes in the surface BZ (k GM/ -mesh). The 
modified offset-I" method designed to treat anisotropic 
systems was employed to perform fc-integration of sur- 
face BZ in the GW part of the code (the method is de- 
scribed in Ref.[H, following the Eq. (53)). The GW 
part of the code, where the self-energy is calculated given 
the eigenfunctions and eigenvalues of the H° generated 
by the LMTO part, is significantly more computation- 
ally demanding than the LMTO calculation. In prac- 
tice, it is computationally prohibitive to calculate the 
self-energy on the same fine k LMT -mesh required for 
the LMTO part. Thus, a rather sophisticated procedure 
that includes several transformations of the self-energy 
S^ ra , between different basis sets has been developed to 
interpolate the self-energy calculated by the GW part of 
the program on a coarse k GM/ -mesh to finer k iMT -mesh 
used by the LMTO part of the program^ 

However, matrix elements of the self-energy (fTJ) be- 
tween states with high energy (> 2 Ry) often cannot be 
interpolated with sufficient accuracy from the h GW -mesh. 
to the k iMTO -mesh2£. Note that a fine k iMTO -mesh is 
required when describing metals. The latter is in part due 
to the long range of the LMTO basis set (e.g. the small- 
est eigenvalue of the overlap matrix can be of the order of 
10 -10 ). In order to overcome this k-interpolation prob- 
lem, the high-energy part (e^ A ,e^? A > E xccut ) of the 
difference between the self-energy and LDA exchange- 
correlation potential, 



av xc — y~ ~ 



c,LDA 
rim 



(3) 



was substituted with a diagonal matrix with the diag- 
onal elements given by linear function of the LDA en- 
ergy, AVjj~ = a + b x e^P A . Here "~" over the sub- 
script denotes that the function is represented in the ba- 
sis of eigenfunctions ip^ A of the LDA Hamiltonian (the 
LDA basis) with eigenvalues s\^? A . The energy cutoff 
parameter E xccut is typically of the order of 2-3 Ry. The 
constants a and b are fitted from calculated AV?? at 
lower energies. The results for the calculated quasipar- 
ticle (QP) energies usually depend weakly on the cutoff 
parameter E xccut or constants a and b. More details on 
the k-interpolation procedure and the method used to 
control its accuracy for bulk calculations could be found 
in Ref. [Hj], section II, subsection G. 

In the LDA basis the optimum Hamiltonian _ff°(k) ([2]) 
reads 

H^(k) = <VC A |ff>k£ A > = ef£ A 8** + AV^(k) . 

(4) 

After modifications of the matrix AVjfj? as outlined 
above, the Hamiltonian H^(\s.) in a form of Eq. @ 



is used in the LMTO part of the program to obtain the 
QP wave functions and energies. 

In the present work, we found that the QSGW method 
requires additional modification when the system ex- 
hibits extended regions with small electronic density, for 
example when modeling a metal/vacuum interface. Ex- 
plicitly, the matrix elements AV?^ in the LDA Hamil- 
tonian of Eq. (T5]) leads to slightly improper mixing be- 
tween occupied LDA states with energies s^? A < Ef 
that are spatially concentrated in the metal region and 
"vacuum" LDA states with energies e£ T f A > ef + $ ex- 
tending to the entire volume of the system (here ep is the 
Fermi energy). As a result, after diagonalization of the 
Hamiltonian (H|), the QP occupied states have small tails 
that decay unphysically slow as a function of distance 
from the metal surface. Several reasons may contribute 
to this slow, non-exponential decay of occupied QP wave 
functions into vacuum, for example: remaining errors in 
the k-interpolation of the AT^j? , non-completeness of the 
LDA basis, numerical errors, etc. 

Development of a general procedure to construct opti- 
mum QSGW Hamiltonian for systems with vacuum re- 
gions in a way that guarantees correct exponential decay 
of occupied QP wave functions in vacuum is beyond the 
scope of this paper. However, in application to the spe- 
cific case of Al, we can overcome this problem by straight- 
forward truncation of the unphysical non-diagonal matrix 
elements AV^ with sfc? A < e F and s£P A > e F + $, 
using the fact that the wave functions of bulk Al are 
rather well described by LDA (see Fig. [T]). Specifically, 
we truncate all non-diagonal matrix elements of AV^ if 
the energies e^n an d £ vK A satisfy following conditions: 



AVft^(k) = AV^~(k) = if 



LDA 



e kfi 



LDA . t? 
£ krh > 



and 



LDA 



> ef 



(5) 



Here E c is a cutoff parameter. The idea of the method 
is to allow the occupied LDA states rh to be mixed by 
matrix AVJfS only with unoccupied LDA states h with 
energy less then sf+E c . Such procedure will prevent oc- 
cupied states from mixing with extended vacuum states 
if E c < $ [at least in the first order of the perturbation 
theory, if consider AVJ?|j in Eq. (|1]) as a perturbation 
to the LDA Hamiltonian]. Note that second condition 
in ([5]) always allows the occupied states to be mixed be- 
tween themselves. In the limit E c —> oo there is no mod- 
ification of the matrix AV^ and the method reduces 
to the standard QSGW approach. In the opposite limit, 
E c — > +0, the method becomes an "unoccupied states 
eigenvalue-only" self-consistent GW method. The "un- 
occupied states eigenvalue-only" means that unoccupied 
subblock of the matrix is diagonal, so the unoc- 

cupied QP wave functions are always equal to the LDA 
wave functions and only QP energies are modified due 
to the diagonal matrix elements AV??. Thus, parameter 
< E c < oo smoothly interpolates between these two 
methods. 
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FIG. 1: (Color online) Density of states (DOS) of bulk Al 
calculated in LDA (dot-dashed line), full QSGW that corre- 
sponds to the limit E c — > oo (solid line), modified QSGW 
with small cutoff parameter E c — 0.1 eV (dotted line), and 
the standard "eigenvalue only" self-consistent GW (dashed 
line). The DOS curves of the full QSGW (E c oo) and 
modified QSGW with E c — 0.1 eV are indistinguishable on 
the figure. 



III. RESULTS AND DISCUSSION 

We begin with analyzing different characteristics of 
bulk Al and Al surfaces as function of parameter E c . 
Figure Q] shows the density of states (DOS) of the bulk 
Al calculated by LDA, full QSGW (that corresponds to 
the limit E c — > oo), modified QSGW as specified by Eq. 
([5]) with small cutoff parameter E c — 0.1 eV, and the 
standard "eigenvalue-only" self-consistent GW methods. 
(In the "eigenvalue-only" self-consistent GW method all 
non-diagonal elements of the GW addition to the LDA 
Hamiltonian are neglected AV?^ = AV? r ~5nm, so the QP 
wave functions are always equal to the LDA ones.) It is 
evident that all three modifications of the GW method 
produce very similar DOS, somewhat different from the 
LDA result. The "eigenvalue-only" DOS is very close to 
the full QSGW DOS with minor deviations in the en- 
ergy range from Ef — 1 eVtoe^ + 3 eV. More impor- 
tantly, the DOS obtained in the full QSGW (solid line) 
and DOS obtained by modified QSGW with small cutoff 
E c = 0.1 eV (dotted line) are almost indistinguishable 
on the figure. This means that the truncation of the 
matrix elements AV^ with arbitrary cutoff parameter 
E c > 0.1 eV practically does not change the bulk Al elec- 
tronic structure. This is an important result suggesting 
that we can safely neglect erroneous non-diagonal ele- 
ments of AV r ^ with e££ A < e F and e£P A > e F + $ in 
surface calculations. 

Next, we turn to the Al/vacuum interface systems. 
Figure shows the electron density averaged over the 
x and y directions for an Al(lll) surface using five dif- 
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FIG. 2: (Color online) The electron density averaged over 
the xy plane (in atomic units, note logarithmic scale) as a 
function of the distance from the center of the vacuum slab, 
z, (in A) for Al(lll)/vacuum interface calculated by modified 
QSGW method using five different values of the E c parameter. 
The number of Al and vacuum layers are Na = 4 and Nv = 
8. The density curves with three smallest values of E c are 
indistinguishable on the figure. 



ferent values of E c . The calculations were performed for 
4 Al and 8 vacuum layers (Na = 4 and Ny = 8). The ex- 
pected well behaved exponential decrease of electron den- 
sity away from the metal surface is seen for E c — 1.36 eV, 
2.72 eV, and 4.08 eV. Importantly, the densities obtained 
using these three E c are indistinguishable. On the other 
hand, for E c above 4.08 eV, the density begins to devi- 
ate from the normal behavior; it sharply increases near 
the center of the vacuum. The density calculated with 
E c = 6.8 eV even increases, at some z, when the distance 
from the metal increases. Similar results, independence 
on E c and correct exponential behavior of density in vac- 
uum for E c < 4 eV, and unphysical behavior for E c > 4 
eV, is found for Al(100) and Al(110) surfaces. 

We note that the unphysical behavior of the electron 
density only occur for small absolute density values, 4-5 
orders of magnitude smaller then the density in metal re- 
gion. Also, the QP energy bands depend only weakly on 
E c : even the energy bands calculated with parameters 
E c = 2.72 eV and 6.80 eV (not shown) almost coincide 
with each other for states with energies less then ep + 2 
eV and begin to deviate slowly for higher energies. In- 
creasing E c from 2.72 eV to 6.8 eV results in an up-shift 
of the QP bands with energy above Ef + 4 eV by a value 
ranged from to 0.2 eV, depending on particular band. 

Figure [3] shows the calculated work function, $, for 
Al(lll) as function of the cutoff parameter E c . The 
value of the work function does not depend on E Cl up 
to E c ~ 4 eV but changes above this threshold. Similar- 
behavior - independence of the work function on E c for 
E c < 4 eV and rapid change above E c ~ 4 eV thresh- 
old, was obtained for Al(100) and Al(110) surfaces. The 
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FIG. 3: The work function calculated for Al(lll) with N A 
4 and Ny — 8 as function of the cutoff parameter E c . 



threshold at E c ~ 4 eV at which point the behavior of 
the work function and electronic density in vacuum both 
sharply changes, is roughly equal to the value of the work 
function E c ~ $ ~ 4.2 — 4.4 eV of Al surfaces, an ad- 
ditional indication that the origin of the error is an im- 
proper mixing of the occupied and vacuum states with 
energies efc£ A > ef + 

Summarizing the results shown in Figures 1-3, we con- 
clude that (1) The DOS of bulk Al docs not depend on 
E c for E c > 0.1 eV; (2) for all three Al surfaces the elec- 
tron density is exponential in vacuum, does not depend 
on E c for E c < 4 eV, and demonstrates unphysical be- 
havior for E c > 4 eV; and (3) for all three Al surfaces 
the value of the work function does not depend on E c for 
E c < 4 eV, and sharply changes for E c > 4 eV. There- 
fore, in the range 0.1 eV < E c < 4 eV, E c is large enough 
for AT^f,£ to include all important matrix elements (at 
least at the level of bulk Al), and simultaneously small 
enough to not include AVJ?? that erroneously mix occu- 
pied LDA states with vacuum LDA states. Importantly, 
the work function in this range does not depend on E c 
and thus could be taken as the true QSGW value of the 
work function. Therefore, in all calculations presented 
below the E c parameter is fixed and set to E c = 2.72 eV. 

The range of applicability of the method described by 
Eq. ([5]) is limited to materials (such as Al) for which 
LDA wave functions are adequate so the matrix elements 
of AV^f with e^P A - e^? A > $ can be neglected. For 
any given metal, this condition can be verified on bulk 
level without performing time consuming surface calcu- 
lations. The method is not applicable to metals (such 
as rf-electron Fe and Cu) for which the matrix elements 
AV~i~j with e^P A - e^ A > $ play significant roles. As 
mentioned above, further efforts are required to develop 
an universal QSGW-derived method applicable to Fe, Cu 
and other metals for which simple truncation of the ma- 
trix elements ([5]) does not work. 

Figure 0] shows the variation of the calculated work 
function, $, as function of slab thickness Na- The calcu- 



lations were performed with LDA and QSGW for unre- 
laxed surfaces using the following parameters: E c = 2.72 
eV, N v = 6, (22 x 22) k iMTO -mesh, and (6 x 6) k GW - 
mesh in the surface BZ (for the more anisotropic Al(110) 
surface we used (22 x 16) k iMTO -mesh, and (6 x 4) k GW - 
mesh). For the LDA calculations we used the Barth- 
Hedin 7 functional. $ oscillates as the slab thickness in- 
creases. These oscillations are well known^ and can be 
attributed to quantum-size effects (QSE). The positions 
of the local maximums of the LDA $ at Na = 5,8, and 11 
for the Al(lll) surface, and minimums at Na = 7 for the 
Al(100) surface and at N A = 5 and 12 for the Al(llO) sur- 
face are in agreement with previous DFT calculations 32 . 
The work functions obtained by the QSGW method show 
similar QSE oscillations. We estimate the uncertainty in 
our calculated for N A = 12 [N A = 14 for Al(110)] values 
of the work function due to the QSE as ±0.03 eV, which 
is larger then uncertainties due to other computational 
parameters like the number of k points in the k GW/ -mesh. 

Delerue et. al^i and Freysoldt et. && found sizable 
rcnormalization of the GW self energy in thin semicon- 
ductor films due to the image potential at the interface; 
this effect is as large as 0.2 eV for the band gap of Si slabs 
with a thickness below 3 nrr*22. This is a finite size effect, 
different from QSE. On the other hand, for metallic films 
the image potential inside the metal slab is well screened, 
so it has only a minor effect on occupied states concen- 
trated within the slab. Since the work function is mostly 
affected by occupied states, we do not expect a signifi- 
cant image potential induced correction to the value of 
the QSGW work function. Furthermore, because there 
is no image potential in the LDA approach, this assump- 
tion is supported by the similar behavior of $ for QSGW 
and LDA as a function of slab thickness, see Fig. 2] 

Results from our work function calculations for three 
Al surfaces are shown in Table I in comparison with ex- 
perimental data and results from other theoretical stud- 
ies. For all three surfaces, our LDA/CA results are rela- 
tively close to those obtained by other groups. The values 
of work function of Al(lll) and Al(100) surfaces calcu- 
lated by different groups using the GGA/BPE method 
also are relatively close to each other. Thus, one can 
conclude that the results obtained using a specific DFT 
functional are converged (within 0.1 eV or better accu- 
racy) for different code implementations. On the other 
hand, Table I shows that the work functions obtained by 
using different DFT functionals could deviate by more 
then 0.1 eV: the LDA/CA values of work functions are 
universally smaller by ~ 0.1 eV then LDA/BH values 
while GGA/PBE values are universally smaller then both 
LDA/BH and LDA/Wagner values by as much as 0.3 eV. 
As mentioned in the introduction, such discrepancies em- 
phasize the need for improved methods if an accuracy of 
0.1 eV or better is required. 

Tabic I shows that the work functions calculated us- 
ing the QSGW method for relaxed Al(lll), Al(100), 
and Al(110) surfaces are equal to 4.17 eV, 4.36 eV, and 
4.19 eV, respectively. We verified that these values do 
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TABLE I: A comparison of experimental values of Al work 
functions (given in eV) for (111), (100), and (110) sur- 
faces with corresponding values calculated by different meth- 
ods: QSGW, G Wo, G W with Al described by jel- 
lium (Go Wo pel]), V IC (GW), LDA/BH, LDA/CA, LDA with 
Wigner interpolation formula^, and GGA/BPE. 



FIG. 4: (Color online) The work function of Al(lll) (top 
panel), Al(100) (middle panel), and Al(110) (bottom panel) 
surfaces calculated by LDA (squares) and QSGW (circles) 
approaches as function of the number of Al layers, Na ■ 



not depend on the particular LDA exchange-correlation 
functional used for the initial iteration by applying both 
LDA/BH and LDA/CA. We found that relaxation of the 
Al surface leads to a less then 0.01 eV shift in the value of 
the QSGW work function: work functions for unrelaxed 
systems are 0.002 eV and 0.007 eV higher for A(lll) and 
Al (100), and 0.008 eV lower for Al (110) surface relative 
to corresponding values for relaxed surfaces. Such small 
effects of surface relaxation agree well with previous DFT 
calculations (see, e.g., Ref. |32| ). 

When compared with data from experimental pho- 
toelectric measurements carried out under ultrahigh 
vacuum^, the work functions obtained using the QSGW 
method for Al(lll), Al(100), and Al(110) surfaces differ 
by 0.07 eV, 0.05 eV, and 0.09 eV, respectively. All three 
differences are less then 0.1 eV and of the order of the sum 
of the theoretical and experimental error bars; we there- 
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fore consider this agreement excellent. Of particular in- 
terest are the differences between different surface faces: 
$(100) - $(111) = 0.19 eV and $(110) - $(111) = 0.02 
eV are both in agreement with the experimental data^ 
$(100) - $(111) = 0.17 eV and $(110) - $(111) = 0.04 
eV. 

Note that all calculated work functions presented in 
Table I (except Go Wo [jel]) follow the increasing trend 
$(111) < $(110) < $(100), in agreement with data. 
This behavior is considered an anomaly; most other 
fee metals instead follow Smoluchowski's rule $(110) < 
$(100) < $(111). The anomaly is caused by an increased 
p-atomic-like character of DOS at the Fermi energy in 
aluminum for the three surfaces, a behavior different 
to that of most fee metals^. We in this context note 
that earlier experiments 4 - reported 0.21 eV and 0.22 eV 
smaller values for Al(lOO) and Al(110) work functions 
compare to that of Ref. [43(. However, Grepstad et al4^ 
suggested that this discrepancy could be due to higher 
impurity concentration, in particular oxygen, in the ear- 
lier experiment. 

Table I shows that Al work functions calculated using 
the G W method differ little (0.01-0.02 eV) from the 
converged QSGW results. The Go Wo results presented 
in Table I correspond to using LDA/BH as starting point 
for GW iterations. Similar 0.01-0.02 eV deviations from 
the converged QSGW results were found for Go Wo when 
instead using LDA/CA. We also note that the conver- 
gence of the GW iterations is fast, meaning that the ini- 
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tial LDA wave functions are close to the converged QP 
wave functions; a conclusion supported by the similarity 
of the QSGW and QSGW(e-only) DOS for bulk Al shown 
in Fig. 1. The substantial differences to previous GoW) 
calculations by Morris et al^i (Go Wo Del] line in Table 
I) and Heinrichsmeier et al 2 ^ (Vz C (GW) line in Table I) 
should therefore neither be attributed to the non-self- 
consistency of the Go Wo method nor to errors associated 
with the choice of particular DFT functional. Instead, 
we propose that the differences are due to the jellium 
approximation employed in both those studie a 24 > 25 . 

For more correlated materials such as Fe or Cu (where 
LDA and GW wave functions overlap less than they do 
in Al) we expect larger deviations of the work functions 
calculated by Go Wo and QSGW methods as well as a 
stronger dependence of the Go Wo results on the DFT 
functional used to calculate Go and Wo. 

IV. SUMMARY 

We have applied the QSGW and Go Wo methods to 
calculate the work functions of Al(lll), Al(100), and 
Al(llO) surfaces. The Go Wo results differ from converged 
QSGW results by less then 0.02 eV and this small differ- 
ence can be attributed to significant overlap of the LDA 
and QP wave functions. The QSGW results are in ex- 
cellent agreement with experimental data taken under 
ultrahigh vacuum conditions. The calculated values of 
the work functions do not depend on the DFT functional 
used for the initial Hamiltonian H°. These results sug- 
gest that QSGW method can be used for reliable and 
accurate calculation of the work functions with accuracy 
of the order of 0.1 eV or better. 



We found that modifications of the original QSGW 
mctho d 10 ' 1:L i 26 are required in order to apply the method 
to the metal/vacuum surface. In particular, special care 
should be taken to control the errors originated form 
(slight) improper mixing of the occupied and vacuum 
states. In some simple cases, such as Al, where LDA 
wave functions are already a good approximation to the 
QP wave functions, simple truncation of corresponding 
matrix elements [see Eq. ((SJ] are enough to control these 
errors. 

The truncation method is not applicable to metals such 
as d-electron Fe and Cu for which the matrix elements 

^■ V Arh witn e kf A ~ £ kfh, A > * P la y significant roles. 
For any given metal, this condition can be verified by 
studying the bulk electronic structure, thus without per- 
forming time consuming surface calculations. Further ef- 
forts are required to develop an universal QSGW-derived 
method applicable to Fe, Cu and other metals for which 
simple truncation of the matrix elements is not adequate. 
V. ACKNOWLEDGEMENT 



We thank Mark van Schilfgaarde for helpful discus- 
sions. This work was supported by the Science of Ex- 
treme Environments LDRD Investment Area at Sandia 
National Laboratories. Sandia is a multiprogram labora- 
tory operated by Sandia Corporation, a Lockheed Mar- 
tin Company, for the United States Department of En- 
ergy under contract DE-AC04-94-AL85000. O.M. and 
S.F. acknowledge the CNMS User support by Oak Ridge 
National Laboratory Division of Scientific User facilities, 
Office of Basic Energy Sciences, U.S. Department of En- 
ergy. 



* Electronic address: sfaleev@mint.ua.edu 

1 D.R.Jennison, P.A. Schultz, D.B. King, and K.R. Zavadil, 
Surface Science 549, 115 (2004). 

2 P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964). 

3 W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965). 

4 H. Kawano, Progress in Sirf. Sci. 83, 1 (2008). 

5 J.P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 
77, 3865 (1996). 

6 D. M. Ceperley and B. J. Alder, Phys. Rev. Lett. 45, 566 
(1980). 

7 U. Von Barth and L. Hedin, J. Phys. C 5, 1629 (1972). 

8 L. Hedin and S. Lundqvist, Solid State Physics (Academic 
Press, New York, 1969), Vol. 23. 

9 M.S. Hybertsen and S.G.Louie, Phys. Rev. Lett. 55, 1418 
(1985). 

M. van Schilfgaarde, T. Kotani, and S. V. Faleev, Phys. 
Rev. Lett. 96, 226402 (2006). 

1 S. V. Faleev, M. van Schilfgaarde, and T. Kotani, Phys. 
Rev. Lett. 93, 126406 (2004). 

2 M.S. Hybertsen and S.G.Louie, Phys. Rev. B 38, 4033 
(1988). 

3 M. Rohlfing, N.-P. Wang, P. Kruger, and J. Pollmann, 



Phys. Rev. Lett. 91,256802 (2003). 

14 I.D. White, R.W. Godby, M.M. Rieger, and R.J. Needs, 
Phys. Rev. Lett. 80, 4265 (1998). 

15 G. Fratesi, G.P Brivio, P. Rinke, and R.W. Godby, Phys. 
Rev. B 68, 195404 (2003). 

16 S. Crampin, Phys. Rev. Lett. 95, 046801 (2005). 

17 P. Rinke, K. Delaney, P. Garcia-Gonzalez, and R.W. 
Godby, Phys. Rev. A 70, 063201 (2004). 

18 J.B. Neaton, M.S. Hybertsen, and S.G.Louie, Phys. Rev. 
Lett. 97, 216405 (2006). 

19 K.S. Thygesen and A. Rubio, Phys. Rev. Lett. 102, 046802 
(2009). 

20 C. Freysoldt, P. Rinke, and M. Schemer, Phys. Rev. Lett. 
103, 056803 (2009). 

21 P. Darancet, A. Ferretti, D. Mayou, and V. Olevano, Phys. 
Rev. B 75, 075102 (2007). 

22 K.S. Thygesen and A. Rubio, I. Chem. Phys. 126, 091101 
(2007). 

23 CD. Spataru, M. S. Hybertsen, S.G. Louie, and A.J. Mol- 
lis, Phys. Rev. B 79, 155110 (2009). 

24 A.J. Morris, M. Stankovski, K. T. Delaney, P. Rinke, 
P. Garcia-Gonzalez, and R.W. Godby, Phys. Rev. B 76, 



8 



25 


155106 (2007). 


35 


E. Wigner, Phys. Rev. 46, 1002 (1934). 


M. Heinrichsmeier, A. Fleszar, W. Hanke, and A. G. 


36 


C.J. Fall, N. Binggeli, and A. Baldereschi, Phys. Rev. B 


26 


Eguiluz, Phys. Rev. B 57, 14974 (1998). 


37 


58, R7544 (1998). 


T. Kotani, M. van Schilfgaarde, and S. V. Faleev, Phys. 


M. Heinrichsmeier, A. Fleszar, and A. G. Eguiluz, Surf. 




Rev. B 76, 165106 (2007). 




Sci. 285, 129 (1993). 


27 


R. W. G. Wyckoff, Crystal Structures, 2nd ed. (Inter- 


38 


J.L.F. Da Silva, C. Stampfl, M. Scheffler, Surf. Sci. 600, 


28 


science, New York, 1963). 


39 


703 (2006). 


M. van Schilfgaarde, T. Kotani, and S. V. Faleev, Phys. 


J. Schochlin, K.P Bohnen, K.M. Ho, Surf. Sci. 324, 113 


29 


Rev. B 74, 245125 (2006). 


40 


(1995). 


C. Frcysoldt, P. Eggert, P. Rinke, A. Schindlmayer, And 


C.J. Fall, N. Binggeli, and A. Baldereschi, Phys. Rev. Lett. 


30 


M. Scheffler, Phys. Reb. B 77, 235428 (2008). 




88, 156802 (2002). 


S. Ismail-Beigi, Phys. Rev. B 73, 233103 (2006). 


41 


A. Kiejna, B.I. Lundqvist, Phys. Rev. B 63, 085405 (2001). 


31 


C.A. Rozzi, D. Varsano, A. Marini, E.K.U. Gross, and A. 


42 


S. J. Sferco, P. Blaha, and K. Schwarz, Phys. Rev. B 76, 




Rubio, Phys. Rev. B 73, 205119 (2006). 




075428 (2007). 


32 


J.L.F. Da Silva, Phys. Rev. B 71, 195416 (2005). 


43 


J. K. Grepstad, P. O. Garland, and B. J. Slagsvold, Surf. 


33 


H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 




Sci. 57, 348 (1976). 


34 


(1976). 


44 


R. M. Eastment and C. H. B. Mee, J. Phys. F: Met. Phys. 


C. Delerue, G. Allan, and M. Lannoo, Phys. Rev. Lett. 90, 




3, 1738 (1973). 




076803 (2003). 





